Ocena wydajności modelu

# library(reticulate)
# use_condaenv("/Users/tomasz/opt/anaconda3/bin/python", required = TRUE)

Testowanie istotności parametrów

Test t-Studenta

Dla każdego parametru \(\theta_i\) modelu GARCH przeprowadzamy test istotności, aby sprawdzić, czy jego wartość jest statystycznie różna od zera. W tym celu używa się testu t-Studenta, który sprawdza hipotezę zerową \[ H_0: \theta_i = 0, \] przeciwko \[ H_1: \theta_i \neq 0. \] Statystyką testową jest \[ t_i = \frac{\hat{\theta_i}}{SE(\hat{\theta_i})}, \] gdzie \(\hat{\theta_i}\) to oszacowana wartość parametru, a \(SE(\hat{\theta_i})\) to jego błąd standardowy. Asymptotycznie \(t_i\) ma rozkład normalny, więc wykonujemy sprawdzenie \[ |t_i| > z_{1-\alpha/2}, \] gdzie \(z_{1-\alpha/2}\) to kwantyl rozkładu normalnego standardowego dla poziomu istotności \(\alpha\).

Alternatywnie, możemy patrzeć na poziom krytyczny \(p\)-value, który jest obliczany jako \[ p = 2 \cdot (1 - \Phi(|t_i|)), \] gdzie \(\Phi\) to funkcja dystrybuanty rozkładu normalnego. Jeśli \(p < \alpha\), to odrzucamy hipotezę zerową i uznajemy, że parametr jest istotny statystycznie.

  • W praktyce, zazwyczaj, gdy \(|t_i| > 2\), to możemy uznać, że parametr jest istotny na poziomie \(\alpha = 0.05\).

Testowanie istotności parametrów

Istotność parametrów dla modelu GARCH(1,1)

  • Wczytujemy dane i obliczamy logarytmiczne stopy zwrotu z akcji Google. Wykres logarytmicznych stóp zwrotu przedstawia zmienność cen akcji w czasie.
googl = pd.read_csv('data/aapl.us.txt', index_col=0, skiprows=2)
googl['log_return'] = np.log1p(googl.iloc[:,0].pct_change())
googl = googl.log_return
googl.index = pd.to_datetime(googl.index, format='%Y-%m-%d')
googl.dropna(inplace=True)

googl.plot(title='Log Returns of Google Stock', figsize=(14, 4))

  • Dopasowujemy model.
model = arch_model(googl, p = 1, q = 1, mean = 'constant', vol = 'GARCH', dist = 't')
model_fit = model.fit(disp='off', cov_type = 'robust')

model_fit.summary()
Constant Mean - GARCH Model Results
Dep. Variable: log_return R-squared: 0.000
Mean Model: Constant Mean Adj. R-squared: 0.000
Vol Model: GARCH Log-Likelihood: 11665.0
Distribution: Standardized Student's t AIC: -23320.0
Method: Maximum Likelihood BIC: -23284.8
No. Observations: 8361
Date: Sat, May 03 2025 Df Residuals: 8360
Time: 15:24:26 Df Model: 1
Mean Model
coef std err t P>|t| 95.0% Conf. Int.
mu -5.9876e-03 3.137e-04 -19.088 3.187e-81 [-6.602e-03,-5.373e-03]
Volatility Model
coef std err t P>|t| 95.0% Conf. Int.
omega 8.3866e-03 5.463e-05 153.514 0.000 [8.280e-03,8.494e-03]
alpha[1] 0.8100 8.104e-02 9.996 1.594e-23 [ 0.651, 0.969]
beta[1] 9.7975e-06 5.900e-03 1.661e-03 0.999 [-1.155e-02,1.157e-02]
Distribution
coef std err t P>|t| 95.0% Conf. Int.
nu 125.9364 1.673 75.284 0.000 [1.227e+02,1.292e+02]


Covariance estimator: robust

Testowanie istotności parametrów

Estymacja parametrów

Maximum Likelihood Estimation (MLE)
Robust Standard Errors (Bollerslev-Wooldridge)

Walidacja założeń modelu

Wykres residuów

  • Szukamy wzorców w resztach modelu, które mogą sugerować, że model nie jest odpowiedni. Wykresy residuów powinny być losowe i nie wykazywać żadnych wyraźnych wzorców.
  • Zgodnie z założeniami reszty powinny być z rozkładu normalnego, widać jednak często występujące ekstremalne wartości.

Walidacja założeń modelu

Autokorelacja residuów - ACF

from statsmodels.graphics.tsaplots import plot_acf

fig, ax = plt.subplots(1, 1, figsize=(14, 4))
plot_acf(model_fit.std_resid, lags=40, ax=ax, alpha=0.05)
ax.set_title('ACF of Standardized Residuals')
plt.show()

Walidacja założeń modelu

Test Ljung-Boxa

Test Ljung-Boxa sprawdza, czy reszty modelu są niezależne. Hipotezy testowe są następujące: \[ H_0: \text{reszty są niezależne} \\ H_1: \text{reszty są skorelowane}. \] Robimy wykres \(p\)-value dla różnych opóźnień \(h\).

from statsmodels.stats.diagnostic import acorr_ljungbox


ljung_box = acorr_ljungbox(model_fit.std_resid, return_df=True)
ljung_box['lb_pvalue'].plot(title='Ljung-Box Test p-values', figsize=(14, 4), style='.')
plt.axhline(y=0.05, color='r', linestyle='--')
plt.xlabel('Lag')
plt.ylabel('p-value')
plt.show()

Miary dopasowania modelu

Akaike Information Criterion (AIC)

\[ AIC = -2 \cdot \log(L) + 2 \cdot k \] gdzie \(L\) to funkcja wiarygodności, a \(k\) to liczba parametrów w modelu. Im mniejsza wartość AIC, tym lepsze dopasowanie modelu.

  • Preferuje modele o mniejszej liczbie parametrów, ale nie karze ich zbytnio za złożoność.
  • Jeśli zależy nam na lepszych zdolnościach predykcyjnych, to lepiej użyć AIC.

Bayesian Information Criterion (BIC)

\[ BIC = -2 \cdot \log(L) + k \cdot \log(n) \] gdzie \(L\) to funkcja wiarygodności, \(k\) to liczba parametrów w modelu, a \(n\) to liczba obserwacji. Podobnie jak AIC, im mniejsza wartość BIC, tym lepsze dopasowanie modelu.

  • Preferuje modele o mniejszej liczbie parametrów, ale karze je bardziej za złożoność.
  • Jeśli zależy nam na lepszym dopasowaniu modelu do danych, to lepiej użyć BIC.

Backtesting modelu GARCH(1,1)

Miary

Mean Absolute Error (MAE)

\[ MAE = \frac{1}{n} \sum_{t=1}^{n} |y_t - \hat{y_t}| \]

Root Mean Squared Error (RMSE)

\[ RMSE = \sqrt{\frac{1}{n} \sum_{t=1}^{n} (y_t - \hat{y_t})^2} \]

  • Wykonujemy prognozę o jeden krok do przodu, a następnie porównujemy prognozowane wartości z rzeczywistymi wartościami. Wykonujemy prognozę na podstawie modelu GARCH(1,1) i obliczamy MAE i RMSE.
from sklearn.metrics import mean_absolute_error, mean_squared_error

# Wykonujemy prognozę na 1 krok do przodu
forecast = model_fit.forecast(horizon=1)
predicted = forecast.variance.values[-1, :][0]
actual = googl.iloc[-1]**2
mae = mean_absolute_error([actual], [predicted])
rmse = np.sqrt(mean_squared_error([actual], [predicted]))
print(f'MAE: {mae}')
print(f'RMSE: {rmse}')
MAE: 0.008448250442643305
RMSE: 0.008448250442643305

ARMA-GARCH w pakiecie R

#| echo: FALSE
#| include: FALSE

library(xts)
library(zoo)
library(rugarch)
library(ggplot2)
library(dplyr)
#| echo: FALSE
#| include: FALSE

googl <- read.csv('data/googl.us.txt', header = TRUE, sep = ",")
googl_xts <- xts(googl[, 5], order.by = as.Date(googl[, 1]))
log_returns <- diff(log(googl_xts))
log_returns <- na.omit(log_returns)
ggplot(log_returns, aes(x = index(log_returns), y = log_returns)) +
  geom_line() +
  labs(title = "Log Returns of Google Stock", x = "Date", y = "Log Returns") +
  theme_minimal()
#| echo: TRUE
spec <- ugarchspec(mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
              variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
              distribution.model = "norm")
fit <- ugarchfit(spec, log_returns)
print(fit)
ggplot(data.frame(fitted = fit@fit$fitted.values, time = index(log_returns)), aes(x = time, y = fitted)) +
  geom_line() +
  labs(title = "Fitted GARCH Model", x = "Date", y = "Fitted Values") +
  theme_minimal()
ggplot(data.frame(residuals = fit@fit$residuals, time = index(log_returns)), aes(x = time, y = residuals)) +
  geom_line() +
  labs(title = "Residuals of GARCH Model", x = "Date", y = "Residuals") +
  theme_minimal()